/*
	The MIT License

	Copyright (c) 2010 IFMO/GameDev Studio

	Permission is hereby granted, free of charge, to any person obtaining a copy
	of this software and associated documentation files (the "Software"), to deal
	in the Software without restriction, including without limitation the rights
	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
	copies of the Software, and to permit persons to whom the Software is
	furnished to do so, subject to the following conditions:

	The above copyright notice and this permission notice shall be included in
	all copies or substantial portions of the Software.

	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
	THE SOFTWARE.
*/

#include "math.h"
#include <stdio.h>
	
/*-----------------------------------------------------------------------------
	QUATERNION.CPP
-----------------------------------------------------------------------------*/

const EQuaternion EQuaternion::kIdentity = EQuaternion(0,0,0,1);


bool EQuaternion::IsEqual ( const EQuaternion &other, float eps ) const
{
	if ( fabs(x - other.x) > eps ) {
		return false;
	}
	if ( fabs(y - other.y) > eps ) {
		return false;
	}
	if ( fabs(z - other.z) > eps ) {
		return false;
	}
	if ( fabs(w - other.w) > eps ) {
		return false;
	}
	return true;
}

EQuaternion& EQuaternion::SetIdentity( void )
{
	x = y = z = 0;
	w = 1;
	return *this;
}

EQuaternion EQuaternion::Conjugate( void ) const
{
	return EQuaternion(-x, -y, -z, w);
}

float EQuaternion::Dot( const EQuaternion &q ) const
{
	return	x * q.x
		+	y * q.y
		+	z * q.z
		+	w * q.w ;
}

EQuaternion EQuaternion::Inverse( void ) const
{
	return Conjugate() / LengthSqr();
}

float EQuaternion::Length( void ) const
{
	return sqrt( Dot(*this) );
}

float EQuaternion::LengthSqr( void ) const
{
	return Dot(*this);
}

EQuaternion EQuaternion::Normalize( void ) const
{
	float len = Length();
	if (len!=0) {
		return EQuaternion(	x/len, y/len, z/len, w/len );
	}
	return *this;
}

//
//	QuatSLerp
//	this function works with non-normalized quaternions,
//	and always returns normalized vector.
//
//	iq = (q*sin((1-t)*omega) + q'*sin(t*omega))/sin(omega),
//	where cos(omega) = inner_product(q,q') 
//
EQuaternion EQuaternion::SLerp( const EQuaternion &target, float factor ) const
{
	EQuaternion	from = *this;
	EQuaternion	to	= target;

	const float DELTA = 0.01f;

	float p1[4];
	double omega, cosom, sinom, scale0, scale1;

	//	cosinus :
	cosom = x*target.x + y*target.y + z*target.z + w*target.w;

	if ( cosom < 0.0 ) { 
		cosom = -cosom;
		p1[0] = - target.x;  p1[1] = - target.y;
		p1[2] = - target.z;  p1[3] = - target.w;
	}
	else {
		p1[0] = target.x;    p1[1] = target.y;
		p1[2] = target.z;    p1[3] = target.w;
	}

	//	angle is bigger than delta :
	if ( (1.0 - cosom) > DELTA ) {
		omega = acos(cosom);
		sinom = sin(omega);
		scale0 = sin((1.0 - factor) * omega) / sinom;
		scale1 = sin(factor * omega) / sinom;
	} else {        
		scale0 = 1.0 - factor;
		scale1 = factor;
	}

	EQuaternion	result(
		(float)( scale0 * x + scale1 * p1[0] ),
		(float)( scale0 * y + scale1 * p1[1] ),
		(float)( scale0 * z + scale1 * p1[2] ),
		(float)( scale0 * w + scale1 * p1[3] ) );

	return result;
}

EMatrix EQuaternion::ToMatrix( void ) const
{
	float	wx, wy, wz;
	float	xx, yy, yz;
	float	xy, xz, zz;
	float	x2, y2, z2;

	x2 = x + x;
	y2 = y + y;
	z2 = z + z;

	xx = x * x2;
	xy = x * y2;
	xz = x * z2;

	yy = y * y2;
	yz = y * z2;
	zz = z * z2;

	wx = w * x2;
	wy = w * y2;
	wz = w * z2;

	EMatrix m;

	m(0,0) = 1.0f - (yy + zz);
	m(1,0) = xy - wz;
	m(2,0) = xz + wy;

	m(0,1) = xy + wz;
	m(1,1) = 1.0f - (xx + zz);
	m(2,1) = yz - wx;

	m(0,2) = xz - wy;
	m(1,2) = yz + wx;
	m(2,2) = 1.0f - (xx + yy);


	return m;
}


void EQuaternion::ToAngles( float &yaw, float &pitch, float &roll ) const
{
	EMatrix	M	=	ToMatrix();
	M.ToAngles( yaw, pitch, roll );
}


void EQuaternion::ToAnglesRad( float &yaw, float &pitch, float &roll ) const
{
	EMatrix	M	=	ToMatrix();
	M.ToAngles( yaw, pitch, roll );

	yaw		=	EMath::Rad(yaw);
	pitch	=	EMath::Rad(pitch);
	roll	=	EMath::Rad(roll);
}


EQuaternion EQuaternion::FromMatrix( const EMatrix &xform )
{
	EMatrix	mat	=	xform.Transpose();

	EQuaternion	q;
	float		trace;
	float		s;
	float		t;
	int     	i;
	int			j;
	int			k;

	static int 	next[ 3 ] = { 1, 2, 0 };

	trace = mat(0,0) + mat(1,1) + mat(2,2);

	if ( trace > 0.0f ) {

		t = trace + 1.0f;
		s = 0.5f / sqrtf( t );

		q.v[3] = s * t;
		q.v[0] = ( mat(2,1) - mat(1,2) ) * s;
		q.v[1] = ( mat(0,2) - mat(2,0) ) * s;
		q.v[2] = ( mat(1,0) - mat(0,1) ) * s;

	} else {

		i = 0;
		if ( mat(1,1) > mat(0,0) ) {
			i = 1;
		}
		if ( mat(2,2) > mat(i,i) ) {
			i = 2;
		}
		j = next[i];
		k = next[j];

		t = ( mat(i,i) - ( mat(j,j) + mat(k,k) ) ) + 1.0f;
		s = 0.5f / sqrt( t );

		q.v[i] = s * t;
		q.v[3] = ( mat(k,j) - mat(j,k) ) * s;
		q.v[j] = ( mat(j,i) + mat(i,j) ) * s;
		q.v[k] = ( mat(k,i) + mat(i,k) ) * s;
	}
	return q;
}


EQuaternion EQuaternion::FromAngles( float yaw, float pitch, float roll )
{
	float	rad_yaw		=	EMath::Rad( yaw	  );
	float	rad_pitch	=	EMath::Rad( pitch );
	float	rad_roll	=	EMath::Rad( roll  );
	return FromAnglesRad( rad_yaw, rad_pitch, rad_roll );
}


EQuaternion EQuaternion::FromAnglesRad( float yaw, float pitch, float roll )
{
	//	TODO : optimize QuatFromAnglesRad
	EQuaternion	qrx		=	RotateAroundAxis( (roll),	EVector(1,0,0) );
	EQuaternion	qry		=	RotateAroundAxis( (pitch),	EVector(0,1,0) );
	EQuaternion	qrz		=	RotateAroundAxis( (yaw),	EVector(0,0,1) );
	EQuaternion	q		=	qrz * qry * qrx;

	return q;
}


EQuaternion EQuaternion::RotateAroundAxis( float angle, const EVector &axis )
{
	//	to return normalized quaternion - normalize axis vector :
	EVector	naxis = axis.Normalize();

	float s = sin(angle/2);
	float c = cos(angle/2);

	return EQuaternion( naxis.x*s, naxis.y*s, naxis.z*s, c );
}


EQuaternion EQuaternion::FromString( const char *str )
{
	EQuaternion p(0,0,0,1);
	sscanf(str, "%f%f%f%f", &p.x, &p.y, &p.z, &p.w);
	return p;
}